# load packages
library(cobalt)
library(vtable)
library(lme4)
library(WeightIt)
library(marginaleffects)
library(patchwork)
library(tidyverse)

# load data
AdminData <- read.csv("BSW/PreparedData.csv")

# define set of official party colours for plotting
CustomColoursGER <- c(
  "AfD" = "#009de0",
  "Greens" = "#409a3c",
  "Left" = "#be3075",
  "CDU" = "#151518",
  "SPD" = "#e3000f",
  "FDP" = "#ffd700",
  "Turnout" = "grey60"
)

# select only municipalities
# & observations with complete party competition
Municips <- AdminData %>% 
  filter(Verwaltungstyp == "Gemeinde" & incomplete_comp == 0)

# assess statistical power of the study
pwr::pwr.t.test(
  n = 423, d = 0.2,
  type = "two.sample", alternative = "two.sided"
)


# =========================================
#  Descriptive Statistics
# =========================================

# names for labels
VarLabs <- data.frame(
  var = c(
    "ChangeTurnout", "ChangeAfD", "ChangeGreens",
    "ChangeFDP", "ChangeSPD", "ChangeCDU", "ChangeLINKE",
    "bsw_option", "east_german",
    "Turnout19", "ShareAfD19", "ShareGreens19",
    "ShareFDP19", "ShareSPD19", "ShareCDU19", "ShareLINKE19",
    "Einwohnerdichte", "Wanderungssaldo",
    "Durchschnittsalter", "Ausländerquote",
    "imputed_foreign", "Hochqualifiziertenquote",
    "SGB.II.Quote", "Kernhaushaltsverschuldung",
    "ChangeTurnoutEU", "ChangeAfDEU", "ChangeGreensEU",
    "ChangeFDPEU", "ChangeSPDEU", "ChangeCDUEU", "ChangeLINKEEU",
    "Turnout24", "TurnoutEU24", "diff_turn_county_to_ep"
  ),
  labs = c(
    "Change Turnout", "Change Vote Share AfD",
    "Change Vote Share Greens", "Change Vote Share FDP",
    "Change Vote Share SPD", "Change Vote Share CDU",
    "Change Vote Share LINKE", "BSW Participation",
    "East Germany", "Turnout 2019", "Vote Share AfD 2019",
    "Vote Share Greens 2019", "Vote Share FDP 2019",
    "Vote Share SPD 2019", "Vote Share CDU 2019",
    "Vote Share LINKE 2019", "Population Density",
    "Migration Balance", "Average Age in Municipality",
    "Share of Foreigners", "Imputed Value for Share Foreigners",
    "Share of Highly Qualified Workers", "Share of Welfare Recipients",
    "Core Budget Debt", "Change Turnout European Election",
    "Change AfD European Election", "Change Greens European Election",
    "Change FDP European Election", "Change SPD European Election",
    "Change CDU European Election", "Change LINKE European Election",
    "Turnout 2024", "Turnout 2024 European Election", 
    "Difference Turnout Council vs. Europe"
  )
)

# select relevant variables + prepare data for table
DescData <- AdminData %>%
  select(
    ChangeTurnout, ChangeAfD, ChangeGreens,
    ChangeFDP, ChangeSPD, ChangeCDU, ChangeLINKE,
    bsw_option, east_german,
    Turnout19, ShareAfD19, ShareGreens19,
    ShareFDP19, ShareSPD19, ShareCDU19, ShareLINKE19,
    Einwohnerdichte, Wanderungssaldo,
    Durchschnittsalter, Ausländerquote,
    imputed_foreign, Hochqualifiziertenquote,
    SGB.II.Quote, Kernhaushaltsverschuldung,
    ChangeTurnoutEU, ChangeAfDEU, ChangeGreensEU,
    ChangeFDPEU, ChangeSPDEU, ChangeCDUEU, ChangeLINKEEU,
    Turnout24, TurnoutEU24, diff_turn_county_to_ep
  ) %>%
  mutate(
    east_german = factor(ifelse(east_german == 1,
      "Yes", "No"
    )),
    imputed_foreign = factor(ifelse(imputed_foreign == 1,
      "Yes", "No"
    ))
  ) 
  
st(DescData,
  summ = c(
    "notNA(x)", "mean(x)", "sd(x)",
    "min(x)", "max(x)"
  ), out = "latex",
  digits = 1, fixed.digits = T, numformat = NA,
  summ.names = c(
    "N", "Mean", "SD",
    "Min", "Max"
  ), fit.page = NA,
  labels = VarLabs,
  title = "Summary Statistics on Full Sample"
)


# ===============================================
# Covariate Balancing Propensity Score Weighting
# ===============================================


# conduct CBPS weighting
# (resulting weight are from Horvitz-Thompson estimator)
CBPSWeighting <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = Municips,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)

# add propensity score + weights to data
Municips <- Municips %>%
  mutate(
    prop_score = CBPSWeighting$ps,
    cbps_weights = CBPSWeighting$weights
  )


# ==================================
# Assumption Check Covariate Balance
# ==================================

# create love plot of balance for means
# (add var names)
Names4Love <- data.frame(
  old_name = c(
    "prop.score",
    "Einwohnerdichte", "Wanderungssaldo",
    "Durchschnittsalter", "Ausländerquote",
    "imputed_foreign", "Hochqualifiziertenquote",
    "SGB.II.Quote", "Kernhaushaltsverschuldung",
    "Turnout19", "ShareLINKE19", "ShareOthers19",
    "east_german"
  ),
  new_name = c(
    "Propensity Score",
    "Population Density", "Migration Balance",
    "Average Age in Municipality", "Share of Foreigners",
    "Imputed Value Foreigners", "Share of Highly Qualified Workers",
    "Share of Welfare Recipients", "Core Budget Debt",
    "Turnout 2019", "Vote Share LINKE 2019",
    "Combined Vote Share Other Parties 2019", "East Germany"
  )
)

# create plot
love.plot(CBPSWeighting,
  stars = "raw", abs = F,
  var.order = "unadjusted",
  drop.distance = F,
  var.names = Names4Love
) +
  geom_vline(xintercept = 0.1, linetype = "dashed") +
  geom_vline(xintercept = -0.1, linetype = "dashed") +
  ggtitle("") +
  scale_color_manual(values = c("grey1", "firebrick3")) +
  xlim(c(-0.8, 0.8)) +
  xlab("Standardized Mean Differences") +
  ylab("Covariate") +
  theme_bw()


# create distribution plot for propensity score
bal.plot(CBPSWeighting,
  var.name = "prop.score",
  which = "both", type = "histogram", bins = 20,
  mirror = T
) +
  ggtitle("") +
  scale_y_continuous(
    breaks = c(-0.10, -0.05, 0, 0.05, 0.10),
    labels = c("0.10", "0.05", "0.00", "0.05", "0.10")
  ) +
  labs(fill = "BSW Participation") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  xlab("Propensity Score") +
  theme_bw()



# create distribution plot for covariates

# Einwohnerdichte
PopDensPlot <- bal.plot(CBPSWeighting,
  var.name = "Einwohnerdichte",
  which = "both", type = "density"
) +
  ggtitle("Population Density") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Density (Total Population / Area in ha)") +
  theme_bw()

# Wanderungssaldo
MigBalPlot <- bal.plot(CBPSWeighting,
  var.name = "Wanderungssaldo",
  which = "both", type = "density"
) +
  ggtitle("Migration Balance") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Migration Balance (per 1,000 Inhabitants)") +
  theme_bw()

# Durchschnittsalter
AverageAgePlot <- bal.plot(CBPSWeighting,
  var.name = "Durchschnittsalter",
  which = "both", type = "density"
) +
  ggtitle("Average Age in Municipality") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Average Age (Years)") +
  theme_bw()

# Ausländerquote
ShareForeignPlot <- bal.plot(CBPSWeighting,
  var.name = "Ausländerquote",
  which = "both", type = "density"
) +
  ggtitle("Share of Foreigners") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Share of Foreigners (%)") +
  theme_bw()

# Imputed Value for Share of Foreigners
ImpForeignPlot <- bal.plot(CBPSWeighting,
  var.name = "imputed_foreign",
  which = "both", type = "density"
) +
  ggtitle("Imputed Values from County Level") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Imputed Value (Share of Foreigners)") +
  theme_bw() 

# Hochqualifiziertenquote
HighQualWorkPlot <- bal.plot(CBPSWeighting,
  var.name = "Hochqualifiziertenquote",
  which = "both", type = "density"
) +
  ggtitle("Share of Highly Qualified Workers") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Employed Residents with Academic Degree (%)") +
  theme_bw()

# SGB.II.Quote
WelfareRecipPlot <- bal.plot(CBPSWeighting,
  var.name = "SGB.II.Quote",
  which = "both", type = "density", alpha = 1
) +
  ggtitle("Share of Welfare Recipients (SGB II)") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Share of Welfare Recipients (%)") +
  theme_bw()

# Kernhaushaltsverschuldung
CoreBudgDebtPlot <- bal.plot(CBPSWeighting,
  var.name = "Kernhaushaltsverschuldung",
  which = "both", type = "density"
) +
  ggtitle("Core Budget Debt of Municipality") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  scale_x_continuous(labels = function(x) x / 1000) +
  labs(fill = "BSW Participation") +
  xlab("Budget Debt (1,000€ per Inhabitant)") +
  theme_bw()

# Turnout 2019
Turnout19Plot <- bal.plot(CBPSWeighting,
  var.name = "Turnout19",
  which = "both", type = "density"
) +
  ggtitle("Turnout in 2019") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Turnout (%)") +
  theme_bw()

# Share LINKE in 2019
ShareLINKEPlot <- bal.plot(CBPSWeighting,
  var.name = "ShareLINKE19",
  which = "both", type = "density"
) +
  ggtitle("Vote Share of the Left in 2019") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Vote Share Left Party (%)") +
  theme_bw()

# Share Other Parties in 2019
ShareOtherPlot <- bal.plot(CBPSWeighting,
  var.name = "ShareOthers19",
  which = "both", type = "density"
) +
  ggtitle("Vote Share of Other Parties in 2019") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("Combined Vote Share of Other Parties (%)") +
  theme_bw()

# East Germany
EastGermanPlot <- bal.plot(CBPSWeighting,
  var.name = "east_german",
  which = "both", type = "density"
) +
  ggtitle("Geographic Location of Municipality") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  xlab("East Germany") +
  theme_bw()

# combine the distribution plots together
ImpForeignPlot +
  EastGermanPlot +
  AverageAgePlot +
  ShareLINKEPlot +
  PopDensPlot +
  WelfareRecipPlot +
  HighQualWorkPlot +
  Turnout19Plot +
  ShareOtherPlot +
  MigBalPlot +
  CoreBudgDebtPlot +
  ShareForeignPlot +
  plot_layout(
    ncol = 3, guides = "collect",
    axis = "collect"
  ) & theme(legend.position = "bottom")

# remove plot objects from environment
rm(
  PopDensPlot, MigBalPlot, AverageAgePlot,
  ShareForeignPlot, ImpForeignPlot, HighQualWorkPlot,
  WelfareRecipPlot, CoreBudgDebtPlot, Turnout19Plot,
  ShareLINKEPlot, ShareOtherPlot, EastGermanPlot
)


# ========================
# Calculate Main Model
# ========================

# CDU
CDUMain <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
CDUMainATE <- avg_comparisons(CDUMain, variables = "bsw_option")

# SPD
SPDMain <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDMainATE <- avg_comparisons(SPDMain, variables = "bsw_option")

# Greens
GreensMain <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensMainATE <- avg_comparisons(GreensMain, variables = "bsw_option")

# FDP
FDPMain <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPMainATE <- avg_comparisons(FDPMain, variables = "bsw_option")

# LINKE
LINKEMain <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKEMainATE <- avg_comparisons(LINKEMain, variables = "bsw_option")

# AfD
AfDMain <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDMainATE <- avg_comparisons(AfDMain, variables = "bsw_option")

# turnout
TurnoutMain <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutMainATE <- avg_comparisons(TurnoutMain, variables = "bsw_option")

# add identifier and bind together
CDUMainATE$party <- "CDU"
SPDMainATE$party <- "SPD"
GreensMainATE$party <- "Greens"
FDPMainATE$party <- "FDP"
LINKEMainATE$party <- "Left"
AfDMainATE$party <- "AfD"
TurnoutMainATE$party <- "Turnout"

MainPlotData <- bind_rows(
  CDUMainATE, SPDMainATE, GreensMainATE,
  FDPMainATE, LINKEMainATE, AfDMainATE,
  TurnoutMainATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUMainATE,
  SPDMainATE,
  GreensMainATE,
  FDPMainATE,
  LINKEMainATE,
  AfDMainATE,
  TurnoutMainATE
)

# create plot
MainATEPlot <- ggplot(MainPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")


# ========================
# Raw Effects Model
# ========================
# Check how much strategic considerations of selection
# into treatment confound results

# CDU
CDUNoWeigh <- lm_weightit(
  ChangeCDU ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
CDUNoWeighATE <- avg_comparisons(CDUNoWeigh, variables = "bsw_option")

# SPD
SPDNoWeigh <- lm_weightit(
  ChangeSPD ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
SPDNoWeighATE <- avg_comparisons(SPDNoWeigh, variables = "bsw_option")

# Greens
GreensNoWeigh <- lm_weightit(
  ChangeGreens ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
GreensNoWeighATE <- avg_comparisons(GreensNoWeigh, variables = "bsw_option")

# FDP
FDPNoWeigh <- lm_weightit(
  ChangeFDP ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
FDPNoWeighATE <- avg_comparisons(FDPNoWeigh, variables = "bsw_option")

# LINKE
LINKENoWeigh <- lm_weightit(
  ChangeLINKE ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
LINKENoWeighATE <- avg_comparisons(LINKENoWeigh, variables = "bsw_option")

# AfD
AfDNoWeigh <- lm_weightit(
  ChangeAfD ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
AfDNoWeighATE <- avg_comparisons(AfDNoWeigh, variables = "bsw_option")

# turnout
TurnoutNoWeigh <- lm_weightit(
  ChangeTurnout ~ bsw_option,
  data = Municips, weightit = NULL,
  cluster = ~Landkreis,
  vcov = "HC0"
)
TurnoutNoWeighATE <- avg_comparisons(TurnoutNoWeigh, variables = "bsw_option")

# add identifier and bind together
CDUNoWeighATE$party <- "CDU"
SPDNoWeighATE$party <- "SPD"
GreensNoWeighATE$party <- "Greens"
FDPNoWeighATE$party <- "FDP"
LINKENoWeighATE$party <- "Left"
AfDNoWeighATE$party <- "AfD"
TurnoutNoWeighATE$party <- "Turnout"

NoWeighPlotData <- bind_rows(
  CDUNoWeighATE, SPDNoWeighATE, GreensNoWeighATE,
  FDPNoWeighATE, LINKENoWeighATE, AfDNoWeighATE,
  TurnoutNoWeighATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUNoWeighATE,
  SPDNoWeighATE,
  GreensNoWeighATE,
  FDPNoWeighATE,
  LINKENoWeighATE,
  AfDNoWeighATE,
  TurnoutNoWeighATE
)

# create plot
NoWeighPlot <- ggplot(NoWeighPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")


# combine the two plots
MainATEPlot +
  ggtitle("Strategic Main Model (CBPS Weights)") +
  xlim(c(-4.5, 2.5)) +
  NoWeighPlot +
  ggtitle("Idiosyncratic Raw Effects Model") +
  xlim(c(-4.5, 2.5)) +
  plot_layout(ncol = 2, guides = "collect", axis = "collect") 
  
rm(MainATEPlot, MainPlotData, NoWeighPlotData, NoWeighPlot,
   AfDMain, AfDNoWeigh, CDUMain, CDUNoWeigh, FDPMain, FDPNoWeigh,
   GreensMain, GreensNoWeigh, LINKEMain, LINKENoWeigh, SPDMain,
   SPDNoWeigh, TurnoutMain, TurnoutNoWeigh)

# ==================
# Robustness Checks
# ==================

# 1) Fit Hierarchical Model with Random Intercepts for Counties
# instead of Clustered Standard Errors

# CDU
CDUHierarch <- lmer(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
CDUHierarchATE <- avg_comparisons(CDUHierarch, variables = "bsw_option", re.form = NA)

# SPD
SPDHierarch <- lmer(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
SPDHierarchATE <- avg_comparisons(SPDHierarch, variables = "bsw_option", re.form = NA)

# Greens
GreensHierarch <- lmer(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
GreensHierarchATE <- avg_comparisons(GreensHierarch, variables = "bsw_option", re.form = NA)

# FDP
FDPHierarch <- lmer(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
FDPHierarchATE <- avg_comparisons(FDPHierarch, variables = "bsw_option", re.form = NA)

# LINKE
LINKEHierarch <- lmer(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
LINKEHierarchATE <- avg_comparisons(LINKEHierarch, variables = "bsw_option", re.form = NA)

# AfD
AfDHierarch <- lmer(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19) + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
AfDHierarchATE <- avg_comparisons(AfDHierarch, variables = "bsw_option", re.form = NA)

# Turnout
TurnoutHierarch <- lmer(
  ChangeTurnout ~ bsw_option +
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german + (1 | Landkreis),
  data = Municips, weights = cbps_weights
)
TurnoutHierarchATE <- avg_comparisons(TurnoutHierarch, variables = "bsw_option", re.form = NA)

# add identifier and bind together
CDUHierarchATE$party <- "CDU"
SPDHierarchATE$party <- "SPD"
GreensHierarchATE$party <- "Greens"
FDPHierarchATE$party <- "FDP"
LINKEHierarchATE$party <- "Left"
AfDHierarchATE$party <- "AfD"
TurnoutHierarchATE$party <- "Turnout"

HierarchPlotData <- bind_rows(
  CDUHierarchATE, SPDHierarchATE, GreensHierarchATE,
  FDPHierarchATE, LINKEHierarchATE, AfDHierarchATE,
  TurnoutHierarchATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))

rm(
  CDUHierarch, CDUHierarchATE,
  SPDHierarch, SPDHierarchATE,
  GreensHierarch, GreensHierarchATE,
  FDPHierarch, FDPHierarchATE,
  LINKEHierarch, LINKEHierarchATE,
  AfDHierarch, AfDHierarchATE,
  TurnoutHierarch, TurnoutHierarchATE
)

# plot ATEs
ggplot(HierarchPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(HierarchPlotData)


# 2) Fit Model Including Observations with Incomplete Competition

# keep observations with incomplete party competition
IncompMunicips <- AdminData %>% 
  filter(Verwaltungstyp == "Gemeinde")

# conduct CBPS weighting
CBPSWeightingIncomp <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = IncompMunicips,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)

# CDU
CDUIncomp <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
CDUIncompATE <- avg_comparisons(CDUIncomp, variables = "bsw_option")

# SPD
SPDIncomp <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDIncompATE <- avg_comparisons(SPDIncomp, variables = "bsw_option")

# Greens
GreensIncomp <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensIncompATE <- avg_comparisons(GreensIncomp, variables = "bsw_option")

# FDP
FDPIncomp <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPIncompATE <- avg_comparisons(FDPIncomp, variables = "bsw_option")

# LINKE
LINKEIncomp <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKEIncompATE <- avg_comparisons(LINKEIncomp, variables = "bsw_option")

# AfD
AfDIncomp <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDIncompATE <- avg_comparisons(AfDIncomp, variables = "bsw_option")

# turnout
TurnoutIncomp <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = IncompMunicips, weightit = CBPSWeightingIncomp,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutIncompATE <- avg_comparisons(TurnoutIncomp, variables = "bsw_option")

# add identifier and bind together
CDUIncompATE$party <- "CDU"
SPDIncompATE$party <- "SPD"
GreensIncompATE$party <- "Greens"
FDPIncompATE$party <- "FDP"
LINKEIncompATE$party <- "Left"
AfDIncompATE$party <- "AfD"
TurnoutIncompATE$party <- "Turnout"

IncompPlotData <- bind_rows(
  CDUIncompATE, SPDIncompATE, GreensIncompATE,
  FDPIncompATE, LINKEIncompATE, AfDIncompATE,
  TurnoutIncompATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUIncomp, CDUIncompATE,
  SPDIncomp, SPDIncompATE,
  GreensIncomp, GreensIncompATE,
  FDPIncomp, FDPIncompATE,
  LINKEIncomp, LINKEIncompATE,
  AfDIncomp, AfDIncompATE,
  TurnoutIncomp, TurnoutIncompATE
)

# plot ATEs
ggplot(IncompPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(IncompPlotData, CBPSWeightingIncomp)


# 3) Fit Model using Entropy Balancing as alternative weighting method

# conduct entropy balancing
EbalWeighting <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = Municips,
  method = "ebal", estimand = "ATE",
  moments = 1, maxit = 10000
)

# CDU
CDUEbal <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
CDUEbalATE <- avg_comparisons(CDUEbal, variables = "bsw_option")

# SPD
SPDEbal <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDEbalATE <- avg_comparisons(SPDEbal, variables = "bsw_option")

# Greens
GreensEbal <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensEbalATE <- avg_comparisons(GreensEbal, variables = "bsw_option")

# FDP
FDPEbal <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPEbalATE <- avg_comparisons(FDPEbal, variables = "bsw_option")

# LINKE
LINKEEbal <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKEEbalATE <- avg_comparisons(LINKEEbal, variables = "bsw_option")

# AfD
AfDEbal <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDEbalATE <- avg_comparisons(AfDEbal, variables = "bsw_option")

# turnout
TurnoutEbal <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = EbalWeighting,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutEbalATE <- avg_comparisons(TurnoutEbal, variables = "bsw_option")

# add identifier and bind together
CDUEbalATE$party <- "CDU"
SPDEbalATE$party <- "SPD"
GreensEbalATE$party <- "Greens"
FDPEbalATE$party <- "FDP"
LINKEEbalATE$party <- "Left"
AfDEbalATE$party <- "AfD"
TurnoutEbalATE$party <- "Turnout"

EbalPlotData <- bind_rows(
  CDUEbalATE, SPDEbalATE, GreensEbalATE,
  FDPEbalATE, LINKEEbalATE, AfDEbalATE,
  TurnoutEbalATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUEbal, CDUEbalATE,
  SPDEbal, SPDEbalATE,
  GreensEbal, GreensEbalATE,
  FDPEbal, FDPEbalATE,
  LINKEEbal, LINKEEbalATE,
  AfDEbal, AfDEbalATE,
  TurnoutEbal, TurnoutEbalATE
)

# plot ATEs
ggplot(EbalPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(EbalPlotData)

# 4) Fit Model that drops observations with imputed values

# filter relevant observations
MunicipsNoImp <- Municips %>%
  filter(imputed_foreign == 0)

# calculate weights
CBPSNoImpWeights <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = MunicipsNoImp,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)

# CDU
CDUNoImp <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  vcov = "asympt"
)
CDUNoImpATE <- avg_comparisons(CDUNoImp, variables = "bsw_option")

# SPD
SPDNoImp <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDNoImpATE <- avg_comparisons(SPDNoImp, variables = "bsw_option")

# Greens
GreensNoImp <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensNoImpATE <- avg_comparisons(GreensNoImp, variables = "bsw_option")

# FDP
FDPNoImp <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPNoImpATE <- avg_comparisons(FDPNoImp, variables = "bsw_option")

# LINKE
LINKENoImp <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKENoImpATE <- avg_comparisons(LINKENoImp, variables = "bsw_option")

# AfD
AfDNoImp <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDNoImpATE <- avg_comparisons(AfDNoImp, variables = "bsw_option")

# turnout
TurnoutNoImp <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MunicipsNoImp, weightit = CBPSNoImpWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutNoImpATE <- avg_comparisons(TurnoutNoImp, variables = "bsw_option")

# add identifier and bind together
CDUNoImpATE$party <- "CDU"
SPDNoImpATE$party <- "SPD"
GreensNoImpATE$party <- "Greens"
FDPNoImpATE$party <- "FDP"
LINKENoImpATE$party <- "Left"
AfDNoImpATE$party <- "AfD"
TurnoutNoImpATE$party <- "Turnout"

NoImpPlotData <- bind_rows(
  CDUNoImpATE, SPDNoImpATE, GreensNoImpATE,
  FDPNoImpATE, LINKENoImpATE, AfDNoImpATE,
  TurnoutNoImpATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUNoImp, CDUNoImpATE,
  SPDNoImp, SPDNoImpATE,
  GreensNoImp, GreensNoImpATE,
  FDPNoImp, FDPNoImpATE,
  LINKENoImp, LINKENoImpATE,
  AfDNoImp, AfDNoImpATE,
  TurnoutNoImp, TurnoutNoImpATE
)

# plot ATEs
ggplot(NoImpPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(NoImpPlotData, MunicipsNoImp, CBPSNoImpWeights)


# 5) Jackknife Saarland because of different voting system
MunicipsNoSaar <- Municips %>%
  filter(Bundesland != "Saarland")

# calculate weights
CBPSNoSaarWeights <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = MunicipsNoSaar,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)

# CDU
CDUNoSaar <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
CDUNoSaarATE <- avg_comparisons(CDUNoSaar, variables = "bsw_option")

# SPD
SPDNoSaar <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDNoSaarATE <- avg_comparisons(SPDNoSaar, variables = "bsw_option")

# Greens
GreensNoSaar <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensNoSaarATE <- avg_comparisons(GreensNoSaar, variables = "bsw_option")

# FDP
FDPNoSaar <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPNoSaarATE <- avg_comparisons(FDPNoSaar, variables = "bsw_option")

# LINKE
LINKENoSaar <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKENoSaarATE <- avg_comparisons(LINKENoSaar, variables = "bsw_option")

# AfD
AfDNoSaar <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDNoSaarATE <- avg_comparisons(AfDNoSaar, variables = "bsw_option")

# turnout
TurnoutNoSaar <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MunicipsNoSaar, weightit = CBPSNoSaarWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutNoSaarATE <- avg_comparisons(TurnoutNoSaar, variables = "bsw_option")

# add identifier and bind together
CDUNoSaarATE$party <- "CDU"
SPDNoSaarATE$party <- "SPD"
GreensNoSaarATE$party <- "Greens"
FDPNoSaarATE$party <- "FDP"
LINKENoSaarATE$party <- "Left"
AfDNoSaarATE$party <- "AfD"
TurnoutNoSaarATE$party <- "Turnout"

NoSaarPlotData <- bind_rows(
  CDUNoSaarATE, SPDNoSaarATE, GreensNoSaarATE,
  FDPNoSaarATE, LINKENoSaarATE, AfDNoSaarATE,
  TurnoutNoSaarATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUNoSaarATE, CDUNoSaar,
  SPDNoSaarATE, SPDNoSaar,
  GreensNoSaarATE, GreensNoSaar,
  FDPNoSaarATE, FDPNoSaar,
  LINKENoSaarATE, LINKENoSaar,
  AfDNoSaarATE, AfDNoSaar,
  TurnoutNoSaarATE, TurnoutNoSaar
)

# plot ATEs
ggplot(NoSaarPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(NoSaarPlotData, MunicipsNoSaar, CBPSNoSaarWeights)

# 6) Sensitivity of Weights
# Outside Common Support of the Propensity Score

# calculate indicator variable for lack of common support
# between treatment and control groups
# using min-max criterion

# calculate common support
Municips %>%
  group_by(bsw_option) %>%
  tally(min(prop_score))

Municips %>%
  group_by(bsw_option) %>%
  tally(max(prop_score))

# set weights to 0 for observations outside common support
Municips <- Municips %>%
  mutate(
    cbps_weights_trimmed = ifelse(
      prop_score < 0.16 | prop_score > 0.84, 0, cbps_weights
    )
  )

# add adapted weights to CBPS object
CBPSWeightsTrimmed <- CBPSWeighting
CBPSWeightsTrimmed$weights <- Municips$cbps_weights_trimmed

# CDU
CDUCommSupp <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis
)
CDUCommSuppATE <- avg_comparisons(CDUCommSupp, variables = "bsw_option")

# SPD
SPDCommSupp <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDCommSuppATE <- avg_comparisons(SPDCommSupp, variables = "bsw_option")

# Greens
GreensCommSupp <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensCommSuppATE <- avg_comparisons(GreensCommSupp, variables = "bsw_option")

# FDP
FDPCommSupp <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPCommSuppATE <- avg_comparisons(FDPCommSupp, variables = "bsw_option")

# LINKE
LINKECommSupp <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKECommSuppATE <- avg_comparisons(LINKECommSupp, variables = "bsw_option")

# AfD
AfDCommSupp <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDCommSuppATE <- avg_comparisons(AfDCommSupp, variables = "bsw_option")

# turnout
TurnoutCommSupp <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSWeightsTrimmed,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutCommSuppATE <- avg_comparisons(TurnoutCommSupp, variables = "bsw_option")

# add identifier and bind together
CDUCommSuppATE$party <- "CDU"
SPDCommSuppATE$party <- "SPD"
GreensCommSuppATE$party <- "Greens"
FDPCommSuppATE$party <- "FDP"
LINKECommSuppATE$party <- "Left"
AfDCommSuppATE$party <- "AfD"
TurnoutCommSuppATE$party <- "Turnout"

CommSuppPlotData <- bind_rows(
  CDUCommSuppATE, SPDCommSuppATE, GreensCommSuppATE,
  FDPCommSuppATE, LINKECommSuppATE, AfDCommSuppATE,
  TurnoutCommSuppATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUCommSuppATE, CDUCommSupp,
  SPDCommSuppATE, SPDCommSupp,
  GreensCommSuppATE, GreensCommSupp,
  FDPCommSuppATE, FDPCommSupp,
  LINKECommSuppATE, LINKECommSupp,
  AfDCommSuppATE, AfDCommSupp,
  TurnoutCommSuppATE, TurnoutCommSupp
)

# plot ATEs
ggplot(CommSuppPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(CommSuppPlotData, CBPSWeightsTrimmed)


# ===================================
# Scope of the Results
# ===================================

# 1) Turnout Effect Without Exactly Concurrent European Election
# (In Thuringia: County Council 26 May 2024 -> European 09 June 2024)

# filter municipalities in Thuringia
# + calculate differences in turnout across elections
Thuringia <- AdminData %>%
  filter(Bundesland == "Thüringen" & incomplete_comp == 0) %>% 
  mutate(diff_turn_county_to_ep = Turnout24 - TurnoutEU24)

# separate control and treatment groups
ThuringiaControl <- Thuringia %>% 
  filter(bsw_option == "No")

ThuringiaTreat <- Thuringia %>% 
  filter(bsw_option == "Yes")

# calculate paired t-tests
t.test(ThuringiaControl$Turnout24, ThuringiaControl$TurnoutEU24,
       var.equal = T, paired = T, alternative = "two.sided")

t.test(ThuringiaTreat$Turnout24, ThuringiaTreat$TurnoutEU24,
       var.equal = T, paired = T, alternative = "two.sided")

# -> Turnout in subsequent EP elections significantly lower in both groups

# compare difference in turnout between the two groups
t.test(ThuringiaTreat$diff_turn_county_to_ep,
       ThuringiaControl$diff_turn_county_to_ep,
       var.equal = F, paired = F, alternative = "greater")

wilcox.test(ThuringiaTreat$diff_turn_county_to_ep,
            ThuringiaControl$diff_turn_county_to_ep,
            paired = F, alternative = "greater")

# generate plot
ThuringiaTurnoutPlot <- ggplot(Thuringia) +
  geom_density(aes(
    x = diff_turn_county_to_ep,
    fill = bsw_option
  ), alpha = 0.5) +
  scale_fill_manual(values = c("#06D6A0", "#FEEFE5")) +
  labs(fill = "BSW Participation") +
  ylab("Density") +
  xlab("Within Unit Change in Turnout (County Council - European)") +
  theme_bw()


# 2) Potential Heterogeneity in East vs. West Germany

# calculate weights balancing East vs. West Germany
CBPSEastWestWeights <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo +
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  by = ~east_german,
  data = Municips,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)


# CDU
CDUEastWest <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
CDUEastWestPred <- avg_comparisons(CDUEastWest, variables = "bsw_option", by = "east_german")
CDUEastWestDiff <- avg_comparisons(CDUEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# SPD
SPDEastWest <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
SPDEastWestPred <- avg_comparisons(SPDEastWest, variables = "bsw_option", by = "east_german")
SPDEastWestDiff <- avg_comparisons(SPDEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# Greens
GreensEastWest <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
GreensEastWestPred <- avg_comparisons(GreensEastWest, variables = "bsw_option", by = "east_german")
GreensEastWestDiff <- avg_comparisons(GreensEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# FDP
FDPEastWest <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
FDPEastWestPred <- avg_comparisons(FDPEastWest, variables = "bsw_option", by = "east_german")
FDPEastWestDiff <- avg_comparisons(FDPEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# LINKE
LINKEEastWest <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
LINKEEastWestPred <- avg_comparisons(LINKEEastWest, variables = "bsw_option", by = "east_german")
LINKEEastWestDiff <- avg_comparisons(LINKEEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# AfD
AfDEastWest <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
AfDEastWestPred <- avg_comparisons(AfDEastWest, variables = "bsw_option", by = "east_german")
AfDEastWestDiff <- avg_comparisons(AfDEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# turnout
TurnoutEastWest <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = Municips, weightit = CBPSEastWestWeights,
  cluster = ~Landkreis
)
TurnoutEastWestPred <- avg_comparisons(TurnoutEastWest, variables = "bsw_option", by = "east_german")
TurnoutEastWestDiff <- avg_comparisons(TurnoutEastWest,
  variables = "bsw_option", by = "east_german",
  hypothesis = "pairwise"
)

# add identifier and bind together
CDUEastWestPred$party <- "CDU"
CDUEastWestDiff$party <- "CDU"
CDUEastWestDiff$contrast <- NA
CDUEastWestDiff$east_german <- NA
SPDEastWestPred$party <- "SPD"
SPDEastWestDiff$party <- "SPD"
SPDEastWestDiff$contrast <- NA
SPDEastWestDiff$east_german <- NA
GreensEastWestPred$party <- "Greens"
GreensEastWestDiff$party <- "Greens"
GreensEastWestDiff$contrast <- NA
GreensEastWestDiff$east_german <- NA
FDPEastWestPred$party <- "FDP"
FDPEastWestDiff$party <- "FDP"
FDPEastWestDiff$contrast <- NA
FDPEastWestDiff$east_german <- NA
LINKEEastWestPred$party <- "Left"
LINKEEastWestDiff$party <- "Left"
LINKEEastWestDiff$contrast <- NA
LINKEEastWestDiff$east_german <- NA
AfDEastWestPred$party <- "AfD"
AfDEastWestDiff$party <- "AfD"
AfDEastWestDiff$contrast <- NA
AfDEastWestDiff$east_german <- NA
TurnoutEastWestPred$party <- "Turnout"
TurnoutEastWestDiff$party <- "Turnout"
TurnoutEastWestDiff$contrast <- NA
TurnoutEastWestDiff$east_german <- NA

PlotDataEastWestPred <- rbind(
  CDUEastWestPred,
  SPDEastWestPred,
  GreensEastWestPred,
  FDPEastWestPred,
  LINKEEastWestPred,
  AfDEastWestPred,
  TurnoutEastWestPred
) %>%
  select(party, term, contrast, east_german, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUEastWestPred, SPDEastWestPred, GreensEastWestPred,
  FDPEastWestPred, LINKEEastWestPred, AfDEastWestPred,
  TurnoutEastWestPred
)

PlotDataEastWestDiff <- rbind(
  CDUEastWestDiff,
  SPDEastWestDiff,
  GreensEastWestDiff,
  FDPEastWestDiff,
  LINKEEastWestDiff,
  AfDEastWestDiff,
  TurnoutEastWestDiff
) %>%
  select(party, term, contrast, east_german, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUEastWestDiff, SPDEastWestDiff, GreensEastWestDiff,
  FDPEastWestDiff, LINKEEastWestDiff, AfDEastWestDiff,
  TurnoutEastWestDiff
)

PlotDataEastWest <- bind_rows(
  PlotDataEastWestPred, PlotDataEastWestDiff
)
rm(PlotDataEastWestPred, PlotDataEastWestDiff)

# adapt factor levels
PlotDataEastWest <- PlotDataEastWest %>%
  mutate(indicator = factor(case_when(
    east_german == 1 ~ "East Germany",
    east_german == 0 ~ "West Germany",
    .default = "First Differences in ATEs (West - East)"
  ), levels = c(
    "West Germany", "East Germany",
    "First Differences in ATEs (West - East)"
  )))

PlotEastWest <- PlotDataEastWest %>%
  filter(indicator == "First Differences in ATEs (West - East)") %>%
  ggplot() +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.2, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("First Differences in ATEs (West - East)") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(
  PlotDataEastWest, CDUEastWest, SPDEastWest, GreensEastWest,
  FDPEastWest, LINKEEastWest, AfDEastWest, TurnoutEastWest, CBPSEastWestWeights
)


# 3) Add county-free cities and county level observations with
# missing municipalities to main dataset (only East Germany)
MixedMunicips <- AdminData %>%
  filter((Verwaltungstyp == "Landkreis" & east_german == 1) |
    Verwaltungstyp == "kreisfreie Stadt")

# adjust data for joining
MixedMunicips$imputed_foreign <- 0
AdjMunicips <- Municips %>%
  select(-cbps_weights, -prop_score, -cbps_weights_trimmed)

# join datasets
MixedMunicips <- bind_rows(
  AdjMunicips, MixedMunicips
)

# remove county Hildburghausen (Greens did not compete)
MixedMunicips <- MixedMunicips %>%
  filter(!(Verwaltungstyp == "Landkreis" & Landkreis == "Hildburghausen, LK"))

# calculate weights
CBPSMixedWeights <- weightit(
  bsw_option ~
    Einwohnerdichte + Wanderungssaldo + 
    Durchschnittsalter + Ausländerquote + imputed_foreign +
    Hochqualifiziertenquote + SGB.II.Quote +
    Kernhaushaltsverschuldung +
    Turnout19 + ShareLINKE19 + ShareOthers19 +
    east_german,
  data = MixedMunicips,
  method = "cbps", estimand = "ATE",
  moments = 1, maxit = 10000
)

# add weights to data
MixedMunicips$cbps_weights <- CBPSMixedWeights$weights


# CDU
CDUMixed <- lm_weightit(
  ChangeCDU ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareCDU19),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
CDUMixedATE <- avg_comparisons(CDUMixed, variables = "bsw_option")

# SPD
SPDMixed <- lm_weightit(
  ChangeSPD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareSPD19),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
SPDMixedATE <- avg_comparisons(SPDMixed, variables = "bsw_option")

# Greens
GreensMixed <- lm_weightit(
  ChangeGreens ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareGreens19),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
GreensMixedATE <- avg_comparisons(GreensMixed, variables = "bsw_option")

# FDP
FDPMixed <- lm_weightit(
  ChangeFDP ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareFDP19),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
FDPMixedATE <- avg_comparisons(FDPMixed, variables = "bsw_option")

# LINKE
LINKEMixed <- lm_weightit(
  ChangeLINKE ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
LINKEMixedATE <- avg_comparisons(LINKEMixed, variables = "bsw_option")

# AfD
AfDMixed <- lm_weightit(
  ChangeAfD ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german +
      ShareAfD19),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
AfDMixedATE <- avg_comparisons(AfDMixed, variables = "bsw_option")

# turnout
TurnoutMixed <- lm_weightit(
  ChangeTurnout ~ bsw_option * (
    Einwohnerdichte + Wanderungssaldo +
      Durchschnittsalter + Ausländerquote + imputed_foreign +
      Hochqualifiziertenquote + SGB.II.Quote +
      Kernhaushaltsverschuldung +
      Turnout19 + ShareLINKE19 + ShareOthers19 +
      east_german),
  data = MixedMunicips, weightit = CBPSMixedWeights,
  cluster = ~Landkreis,
  vcov = "asympt"
)
TurnoutMixedATE <- avg_comparisons(TurnoutMixed, variables = "bsw_option")

# add identifier and bind together
CDUMixedATE$party <- "CDU"
SPDMixedATE$party <- "SPD"
GreensMixedATE$party <- "Greens"
FDPMixedATE$party <- "FDP"
LINKEMixedATE$party <- "Left"
AfDMixedATE$party <- "AfD"
TurnoutMixedATE$party <- "Turnout"

MixedPlotData <- bind_rows(
  CDUMixedATE, SPDMixedATE, GreensMixedATE,
  FDPMixedATE, LINKEMixedATE, AfDMixedATE,
  TurnoutMixedATE
) %>%
  select(party, term, contrast, estimate, conf.low, conf.high) %>%
  mutate(party = factor(party,
                        levels = c("Left", "CDU", "SPD", "FDP",
                                   "Greens", "AfD", "Turnout")))
rm(
  CDUMixedATE, CDUMixed,
  SPDMixedATE, SPDMixed,
  GreensMixedATE, GreensMixed,
  FDPMixedATE, FDPMixed,
  LINKEMixedATE, LINKEMixed,
  AfDMixedATE, AfDMixed,
  TurnoutMixedATE, TurnoutMixed
)

# plot ATEs
MixedPlot <- ggplot(MixedPlotData) +
  geom_errorbar(aes(
    x = estimate, y = party,
    xmin = conf.low, xmax = conf.high,
    colour = party
  ), linewidth = 1.5, width = 0.4) +
  geom_point(aes(
    x = estimate, y = party,
    colour = party
  ), size = 3.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_color_manual(values = CustomColoursGER) +
  xlab("Average Treatment Effect") +
  ylab("Outcome") +
  theme_bw() +
  theme(legend.position = "none")

rm(MixedPlotData, CBPSMixedWeights, AdjMunicips)


# 4) Compare results with those of European Elections

# filter treated municipalities + county-free cities
MunicipsTreat <- MixedMunicips %>%
  filter(bsw_option == "Yes")

# transform data into long format for plotting
MunicipsTreatLong <- MunicipsTreat %>% 
  select(ChangeCDU, ChangeCDUEU, 
         ChangeSPD, ChangeSPDEU,
         ChangeGreens, ChangeGreensEU,
         ChangeLINKE, ChangeLINKEEU,
         ChangeFDP, ChangeFDPEU,
         ChangeAfD, ChangeAfDEU,
         ChangeTurnout, ChangeTurnoutEU) %>%
  pivot_longer(cols = everything(),
               names_to = "election",
               values_to = "change") %>% 
  mutate(party = factor(rep(c("CDU", "CDU", 
                              "SPD", "SPD",
                              "Greens", "Greens",
                              "Left", "Left",
                              "FDP", "FDP",
                              "AfD", "AfD",
                              "Turnout", "Turnout"),
                            227))) %>% 
  relocate(party) %>% 
  mutate(election = factor(rep(c("Local", "EU"), 1589)))

# generate plot
EUChangePlot <- ggplot(MunicipsTreatLong) +
  geom_density(aes(x = change, fill = election), alpha = 0.5) +
  facet_wrap(~party, ncol = 4, scales = "free_y") +
  scale_fill_manual(values = c("grey1", "#FEEFE5")) +
  xlim(c(-25, 25)) +
  labs(fill = "Election Type") +
  ylab("Density") +
  xlab("Change in Outcome (%)") +
  theme_bw()

# combine plots for external validation
Design <- "
112
333
444
"

ThuringiaTurnoutPlot +   
  ggtitle("Distributions of Changes in Turnout Across Non-Concurrent Elections (Thuringia)") +
  guide_area() +
  (PlotEastWest +
  ggtitle("Heterogeneity West and East Germany") |
  MixedPlot +
  ggtitle("Including County-Free Cities and Rural Counties")) +
  EUChangePlot +
  ggtitle("Distributions of Changes in Outcomes in County Council and European Elections") +
  plot_layout(design = Design, guides = "collect")